/******************************/
/** UNEMPLOYED - CSS problem **/
/******************************/



double focU_E(const double xval, void * params){
    
    struct paramsgolden *focparamsU= (struct paramsgolden *) params;
    
    double *expU_UiNE_JNE        = focparamsU->expU_UiNE_JNE_ST;
    double *expU_UnNE_JNE        = focparamsU->expU_UnNE_JNE_ST;
    double *expU_UiE_JE          = focparamsU->expU_UiE_JE_ST;
    double *expU_UnE_JE          = focparamsU->expU_UnE_JE_ST;
    
    double *expU_UiNE            = focparamsU->expU_UiNE_ST;
    double *expU_UnNE            = focparamsU->expU_UnNE_ST;
    double *expU_UiE             = focparamsU->expU_UiE_ST;
    double *expU_UnE             = focparamsU->expU_UnE_ST;
    
    double *expU_WWE             = focparamsU->expU_WWE_ST;
    double *expU_WWNE            = focparamsU->expU_WWNE_ST;
    
    double *expU_WWNE_JNE        = focparamsU->expU_WWNE_JNE_ST;
    double *expU_WWE_JE          = focparamsU->expU_WWE_JE_ST;
    
#if OPT_compute_C_flow == 1
    double *C_expU_UiNE_JNE        = focparamsU->C_expU_UiNE_JNE_ST;
    double *C_expU_UnNE_JNE        = focparamsU->C_expU_UnNE_JNE_ST;
    double *C_expU_UiE_JE          = focparamsU->C_expU_UiE_JE_ST;
    double *C_expU_UnE_JE          = focparamsU->C_expU_UnE_JE_ST;
    
    double *C_expU_UiNE            = focparamsU->C_expU_UiNE_ST;
    double *C_expU_UnNE            = focparamsU->C_expU_UnNE_ST;
    double *C_expU_UiE             = focparamsU->C_expU_UiE_ST;
    double *C_expU_UnE             = focparamsU->C_expU_UnE_ST;
    
    double *C_expU_WWE             = focparamsU->C_expU_WWE_ST;
    double *C_expU_WWNE            = focparamsU->C_expU_WWNE_ST;
    
    double *C_expU_WWNE_JNE        = focparamsU->C_expU_WWNE_JNE_ST;
    double *C_expU_WWE_JE          = focparamsU->C_expU_WWE_JE_ST;
#endif
    
#if SEA == 1
    double *expU_UiE_JEi         = focparamsU->expU_UiE_JEi_ST;
    double *expU_UiNE_JNEi       = focparamsU->expU_UiNE_JNEi_ST;
    double *expU_WWNE_JNEi       = focparamsU->expU_WWNE_JNEi_ST;
    double *expU_WWE_JEi         = focparamsU->expU_WWE_JEi_ST;
#if OPT_compute_C_flow == 1
    double *C_expU_UiE_JEi         = focparamsU->C_expU_UiE_JEi_ST;
    double *C_expU_UiNE_JNEi       = focparamsU->C_expU_UiNE_JNEi_ST;
    double *C_expU_WWNE_JNEi       = focparamsU->C_expU_WWNE_JNEi_ST;
    double *C_expU_WWE_JEi         = focparamsU->C_expU_WWE_JEi_ST;
#endif
#endif
    
    int igridfoc=focparamsU->igridST; int tgridfoc=focparamsU->tgridST;
    int typeU=focparamsU->Utype;
    int excluded=focparamsU->excludedST;
    
    double cons,xgrid,dxgrid,valfoc;
    double tempWWNE, tempWWE, tempUnNE, tempUnE, tempUiNE, tempUiE, tempWWE_JE, tempUiE_JE, tempUiNE_JNE, tempUnE_JE, tempUnNE_JNE, tempWWNE_JNE, tempWWNE_JE, tempJE, tempJNE, tempUiE_JEi, tempUiNE_JNEi,tempJEi, tempJNEi, tempWWNE_JNEi, tempWWE_JEi;
    double C_tempWWNE, C_tempWWE, C_tempUnNE, C_tempUnE, C_tempUiNE, C_tempUiE, C_tempWWE_JE, C_tempUiE_JE, C_tempUiNE_JNE, C_tempUnE_JE, C_tempUnNE_JNE, C_tempWWNE_JNE, C_tempWWNE_JE, C_tempJE, C_tempJNE, C_tempUiE_JEi, C_tempUiNE_JNEi,C_tempJEi, C_tempJNEi, C_tempWWNE_JNEi, C_tempWWE_JEi;
    double testval1, testval2, testval3, testval4;
    double searcheffWW, searcheffJ;
    int ixgrid, iverif, idiffgrid, inxU, inxU1;
    
    
    // CONSUMPTION
    cons   = findconsU(xval, igridfoc, tgridfoc, typeU);
    xgrid  = invexpspace(xval,Gridmin,Gridmax,Echelle1,maxgrid);
    dxgrid = weightinter(xval, xgrid, &ixgrid, grid, maxgrid);
    if ((xval>(Gridmax))){printf("focU: (xval>(Gridmax)) %20.15f\n",xval);getchar();}
    if ((ixgrid>=(maxgrid-1))){printf("focU: ixgrid>=(maxgrid-1) %d\n",ixgrid);getchar();}
    if ((cons<=0.0)){printf("focU: consoFOC<0.0 %20.15f\t%20.15f\n",xval,cons);getchar();}
    if ((xval<(0.0))){printf("focU: (xval<(0.0)) %20.15f\n",xval);getchar();}
    
    
    // CONTINUATION VALUE //
    inxU  = inx(ixgrid,tgridfoc);
    inxU1 = inx((ixgrid+1),tgridfoc);
    
    tempUnNE        = inter1d(dxgrid,expU_UnNE[inxU],expU_UnNE[inxU1]);
    tempUnE         = inter1d(dxgrid,expU_UnE[inxU],expU_UnE[inxU1]);
    tempUiNE        = inter1d(dxgrid,expU_UiNE[inxU],expU_UiNE[inxU1]);
    tempUiE         = inter1d(dxgrid,expU_UiE[inxU],expU_UiE[inxU1]);
    tempWWNE        = inter1d(dxgrid,expU_WWNE[inxU],expU_WWNE[inxU1]);
    tempWWE         = inter1d(dxgrid,expU_WWE[inxU],expU_WWE[inxU1]);
#if OPT_compute_C_flow == 1
    C_tempUnNE        = inter1d(dxgrid,C_expU_UnNE[inxU],C_expU_UnNE[inxU1]);
    C_tempUnE         = inter1d(dxgrid,C_expU_UnE[inxU],C_expU_UnE[inxU1]);
    C_tempUiNE        = inter1d(dxgrid,C_expU_UiNE[inxU],C_expU_UiNE[inxU1]);
    C_tempUiE         = inter1d(dxgrid,C_expU_UiE[inxU],C_expU_UiE[inxU1]);
    C_tempWWNE        = inter1d(dxgrid,C_expU_WWNE[inxU],C_expU_WWNE[inxU1]);
    C_tempWWE         = inter1d(dxgrid,C_expU_WWE[inxU],C_expU_WWE[inxU1]);
#endif
    tempUiNE_JNE    = inter1d(dxgrid,expU_UiNE_JNE[inxU],expU_UiNE_JNE[inxU1]);
    tempUnNE_JNE    = inter1d(dxgrid,expU_UnNE_JNE[inxU],expU_UnNE_JNE[inxU1]);
    tempWWNE_JNE    = inter1d(dxgrid,expU_WWNE_JNE[inxU],expU_WWNE_JNE[inxU1]);
    tempUiE_JE      = inter1d(dxgrid,expU_UiE_JE[inxU],expU_UiE_JE[inxU1]);
    tempUnE_JE      = inter1d(dxgrid,expU_UnE_JE[inxU],expU_UnE_JE[inxU1]);
    tempWWE_JE      = inter1d(dxgrid,expU_WWE_JE[inxU],expU_WWE_JE[inxU1]);
#if OPT_compute_C_flow == 1
    C_tempUiNE_JNE    = inter1d(dxgrid,C_expU_UiNE_JNE[inxU],C_expU_UiNE_JNE[inxU1]);
    C_tempUnNE_JNE    = inter1d(dxgrid,C_expU_UnNE_JNE[inxU],C_expU_UnNE_JNE[inxU1]);
    C_tempWWNE_JNE    = inter1d(dxgrid,C_expU_WWNE_JNE[inxU],C_expU_WWNE_JNE[inxU1]);
    C_tempUiE_JE      = inter1d(dxgrid,C_expU_UiE_JE[inxU],C_expU_UiE_JE[inxU1]);
    C_tempUnE_JE      = inter1d(dxgrid,C_expU_UnE_JE[inxU],C_expU_UnE_JE[inxU1]);
    C_tempWWE_JE      = inter1d(dxgrid,C_expU_WWE_JE[inxU],C_expU_WWE_JE[inxU1]);
#endif
#if SEA == 1
    tempUiE_JEi     = inter1d(dxgrid,expU_UiE_JEi[inxU],expU_UiE_JEi[inxU1]);
    tempUiNE_JNEi   = inter1d(dxgrid,expU_UiNE_JNEi[inxU],expU_UiNE_JNEi[inxU1]);
    tempWWNE_JNEi   = inter1d(dxgrid,expU_WWNE_JNEi[inxU],expU_WWNE_JNEi[inxU1]);
    tempWWE_JEi     = inter1d(dxgrid,expU_WWE_JEi[inxU],expU_WWE_JEi[inxU1]);
#if OPT_compute_C_flow == 1
    C_tempUiE_JEi     = inter1d(dxgrid,C_expU_UiE_JEi[inxU],C_expU_UiE_JEi[inxU1]);
    C_tempUiNE_JNEi   = inter1d(dxgrid,C_expU_UiNE_JNEi[inxU],C_expU_UiNE_JNEi[inxU1]);
    C_tempWWNE_JNEi   = inter1d(dxgrid,C_expU_WWNE_JNEi[inxU],C_expU_WWNE_JNEi[inxU1]);
    C_tempWWE_JEi     = inter1d(dxgrid,C_expU_WWE_JEi[inxU],C_expU_WWE_JEi[inxU1]);
#endif
#endif
    
    
    // STARTING 2D searching effort //
    struct searchfocWwhenE paramsXX;
    
    paramsXX.tempWWNEval = tempWWNE;
    paramsXX.tempUiNEval = tempUiNE;
    paramsXX.tempUnNEval = tempUnNE;
    
    paramsXX.tempWWEval = tempWWE;
    paramsXX.tempUiEval = tempUiE;
    paramsXX.tempUnEval = tempUnE;
    
#if SEA == 1
    paramsXX.tempUiNE_JNEival   = tempUiNE_JNEi;
    paramsXX.tempUiE_JEival     = tempUiE_JEi;
    paramsXX.tempWWNE_JNEival   = tempWWNE_JNEi;
    paramsXX.tempWWE_JEival     = tempWWE_JEi;
#endif
    
    paramsXX.tempWWNE_JNEval = tempWWNE_JNE;
    paramsXX.tempUiNE_JNEval = tempUiNE_JNE;
    paramsXX.tempUnNE_JNEval = tempUnNE_JNE;
    
    paramsXX.tempWWE_JEval = tempWWE_JE;
    paramsXX.tempUiE_JEval = tempUiE_JE;
    paramsXX.tempUnE_JEval = tempUnE_JE;
    
    paramsXX.tgridparXX = tgridfoc;
    paramsXX.igridparXX = igridfoc;
    paramsXX.typeUU     = typeU;
    paramsXX.excludedU  = excluded;
    
    testval1 = searchWwhenEknownfun(0.0, &paramsXX);
    testval2 = searchWwhenEknownfun(TOLZbrent, &paramsXX);
    testval3 = searchWwhenEknownfun(EffortMaxW, &paramsXX);
    testval4 = searchWwhenEknownfun(EffortMaxW-0.01, &paramsXX);
    
    iverif = 0;
    
    if(testval1 > testval2 & testval1 < 0.0) {
        searcheffWW = 0.0;
        searchWwhenEknownfun(searcheffWW, &paramsXX);   // see if you can remove after.
        searcheffJ = paramsXX.searchJoutST;
        iverif = 100;
    }
    if(testval1 < testval2 & testval1 > 0.0) {
        searcheffWW = 0.0;
        searchWwhenEknownfun(searcheffWW, &paramsXX);   // see if you can remove after.
        searcheffJ = paramsXX.searchJoutST;
        iverif = 100;
    }
    
    if(testval3 > testval4 & testval3 < 0.0) {
        searcheffWW = EffortMaxW;
        searchWwhenEknownfun(searcheffWW, &paramsXX);   // see if you can remove after.
        searcheffJ = paramsXX.searchJoutST;
        printf("There is a bug:: EffortmaxW reached \n"); getchar();
        iverif = 100;
    }
    if(testval3 < testval4 & testval3 > 0.0) {
        searcheffWW = EffortMaxW;
        searchWwhenEknownfun(searcheffWW, &paramsXX);   // see if you can remove after.
        searcheffJ = paramsXX.searchJoutST;
        printf("There is a bug:: EffortmaxW reached \n"); getchar();
        iverif = 100;
    }
    
    if(iverif == 0) {
        searcheffWW = zbrentNEW(searchWwhenEknownfun,TOLZbrent,EffortMaxW,TOLZbrent,&paramsXX);
        searchWwhenEknownfun(searcheffWW, &paramsXX);   // see if you can remove after.
        searcheffJ = paramsXX.searchJoutST;
    }
    
    
    // IF SEARCH EFFORT CAN NOT ADJUST //
    if(typeU == 0){
#if holdfixedSE
        if(excluded == 1){searcheffJ = fixedSEUiE[inx(igridfoc, tgridfoc)];}
        if(excluded == 0){searcheffJ = fixedSEUiNE[inx(igridfoc, tgridfoc)];}
#endif
#if holdfixedSW
        if(excluded == 1){searcheffWW = fixedSWUiE[inx(igridfoc, tgridfoc)];}
        if(excluded == 0){searcheffWW = fixedSWUiNE[inx(igridfoc, tgridfoc)];}
#endif
    }
    
    
    focparamsU->searchWoutST=searcheffWW;
    focparamsU->searchJoutST=searcheffJ;
    
    
    
    // COMPUTE THE CONTINUATION VALUE //
    if(excluded == 0) { // not excluded
        if(typeU == 0) {    // insured
#if NOPOLICY == 1
            valfoc = utilc(cons,0) + disutilityW(searcheffWW) + disutilityE(searcheffJ) + betapar*(
                                                                                                   piW(searcheffWW)*((1.0-piE(searcheffJ))*tempWWNE+piE(searcheffJ)*tempWWNE_JNE)+
                                                                                                   (1.0-piW(searcheffWW))*(piE(searcheffJ)*((1.0-pLTpar)*tempUiNE_JNE + pLTpar*tempUnNE_JNE) + (1.0-piE(searcheffJ))*((1.0-pLTpar)*tempUiNE + pLTpar*tempUnNE)));
#if OPT_compute_C_flow == 1
            focparamsU->C_val = utilc(cons,0) + betapar*(
                                                         piW(searcheffWW)*((1.0-piE(searcheffJ))*C_tempWWNE+piE(searcheffJ)*C_tempWWNE_JNE)+
                                                         (1.0-piW(searcheffWW))*(piE(searcheffJ)*((1.0-pLTpar)*C_tempUiNE_JNE + pLTpar*C_tempUnNE_JNE) + (1.0-piE(searcheffJ))*((1.0-pLTpar)*C_tempUiNE + pLTpar*C_tempUnNE)));
#endif
#endif
            
#if SEA == 1
            valfoc = utilc(cons,0) + disutilityW(searcheffWW) + disutilityE(searcheffJ) + betapar*(
                                                                                                   piW(searcheffWW)*((1.0-piE(searcheffJ))*tempWWNE+piE(searcheffJ)*((1.0-pLTpar)*tempWWNE_JNEi + pLTpar*tempWWNE_JNE))+
                                                                                                   (1.0-piW(searcheffWW))*(piE(searcheffJ)*((1.0-pLTpar)*tempUiNE_JNEi + pLTpar*tempUnNE_JNE) + (1.0-piE(searcheffJ))*((1.0-pLTpar)*tempUiNE + pLTpar*tempUnNE)));
#if OPT_compute_C_flow == 1
            focparamsU->C_val = utilc(cons,0) + betapar*(
                                                         piW(searcheffWW)*((1.0-piE(searcheffJ))*C_tempWWNE+piE(searcheffJ)*((1.0-pLTpar)*C_tempWWNE_JNEi + pLTpar*C_tempWWNE_JNE))+
                                                         (1.0-piW(searcheffWW))*(piE(searcheffJ)*((1.0-pLTpar)*C_tempUiNE_JNEi + pLTpar*C_tempUnNE_JNE) + (1.0-piE(searcheffJ))*((1.0-pLTpar)*C_tempUiNE + pLTpar*C_tempUnNE)));
#endif
#endif
        }
        if(typeU == 1) {    // not insured
            valfoc = utilc(cons,0) + disutilityW(searcheffWW) + disutilityE(searcheffJ) + betapar*(
                                                                                                   piW(searcheffWW)*((1.0-piE(searcheffJ))*tempWWNE+piE(searcheffJ)*tempWWNE_JNE)+
                                                                                                   (1.0-piW(searcheffWW))*(piE(searcheffJ)*(tempUnNE_JNE)+(1.0-piE(searcheffJ))*(tempUnNE)));
#if OPT_compute_C_flow == 1
            focparamsU->C_val = utilc(cons,0) + betapar*(
                                                         piW(searcheffWW)*((1.0-piE(searcheffJ))*C_tempWWNE+piE(searcheffJ)*C_tempWWNE_JNE)+
                                                         (1.0-piW(searcheffWW))*(piE(searcheffJ)*(C_tempUnNE_JNE)+(1.0-piE(searcheffJ))*(C_tempUnNE)));
#endif
        }
    }
    
    if(excluded == 1) { // excluded
        if(typeU == 0) {    // insured
#if NOPOLICY == 1
            valfoc = utilc(cons,0)+ disutilityW(searcheffWW) + disutilityE(searcheffJ) + betapar*(
                                                                                                  phipar*(piW(searcheffWW)*((1.0-piE(searcheffJ))*tempWWNE+piE(searcheffJ)*tempWWNE_JNE)+
                                                                                                          (1.0-piW(searcheffWW))*(piE(searcheffJ)*((1.0-pLTpar)*tempUiNE_JNE + pLTpar*tempUnNE_JNE) + (1.0-piE(searcheffJ))*((1.0-pLTpar)*tempUiNE + pLTpar*tempUnNE)))
                                                                                                  +(1.0 - phipar)*(piW(searcheffWW)*((1.0-piE(searcheffJ))*tempWWE+piE(searcheffJ)*tempWWE_JE)+
                                                                                                                   (1.0-piW(searcheffWW))*(piE(searcheffJ)*((1.0-pLTpar)*tempUiE_JE + pLTpar*tempUnE_JE) + (1.0-piE(searcheffJ))*((1.0-pLTpar)*tempUiE + pLTpar*tempUnE))));
#if OPT_compute_C_flow == 1
            focparamsU->C_val = utilc(cons,0) +  betapar*(
                                                          phipar*(piW(searcheffWW)*((1.0-piE(searcheffJ))*C_tempWWNE+piE(searcheffJ)*C_tempWWNE_JNE)+
                                                                  (1.0-piW(searcheffWW))*(piE(searcheffJ)*((1.0-pLTpar)*C_tempUiNE_JNE + pLTpar*C_tempUnNE_JNE) + (1.0-piE(searcheffJ))*((1.0-pLTpar)*C_tempUiNE + pLTpar*C_tempUnNE)))
                                                          +(1.0 - phipar)*(piW(searcheffWW)*((1.0-piE(searcheffJ))*C_tempWWE+piE(searcheffJ)*C_tempWWE_JE)+
                                                                           (1.0-piW(searcheffWW))*(piE(searcheffJ)*((1.0-pLTpar)*C_tempUiE_JE + pLTpar*C_tempUnE_JE) + (1.0-piE(searcheffJ))*((1.0-pLTpar)*C_tempUiE + pLTpar*C_tempUnE))));
#endif
#endif
            
#if SEA == 1
            valfoc = utilc(cons,0)+ disutilityW(searcheffWW) + disutilityE(searcheffJ) + betapar*(
                                                                                                  phipar*(piW(searcheffWW)*((1.0-piE(searcheffJ))*tempWWNE+piE(searcheffJ)*((1.0-pLTpar)*tempWWNE_JNEi + pLTpar*tempWWNE_JNE))+
                                                                                                          (1.0-piW(searcheffWW))*(piE(searcheffJ)*((1.0-pLTpar)*tempUiNE_JNEi + pLTpar*tempUnNE_JNE) + (1.0-piE(searcheffJ))*((1.0-pLTpar)*tempUiNE + pLTpar*tempUnNE)))
                                                                                                  +(1.0-phipar)*(piW(searcheffWW)*((1.0-piE(searcheffJ))*tempWWE+piE(searcheffJ)*((1.0-pLTpar)*tempWWE_JEi + pLTpar*tempWWE_JE))+
                                                                                                                 (1.0-piW(searcheffWW))*(piE(searcheffJ)*((1.0-pLTpar)*tempUiE_JEi + pLTpar*tempUnE_JE) + (1.0-piE(searcheffJ))*((1.0-pLTpar)*tempUiE + pLTpar*tempUnE))));
#if OPT_compute_C_flow == 1
            focparamsU->C_val = utilc(cons,0)+  betapar*(
                                                         phipar*(piW(searcheffWW)*((1.0-piE(searcheffJ))*C_tempWWNE+piE(searcheffJ)*((1.0-pLTpar)*C_tempWWNE_JNEi + pLTpar*C_tempWWNE_JNE))+
                                                                 (1.0-piW(searcheffWW))*(piE(searcheffJ)*((1.0-pLTpar)*C_tempUiNE_JNEi + pLTpar*C_tempUnNE_JNE) + (1.0-piE(searcheffJ))*((1.0-pLTpar)*C_tempUiNE + pLTpar*C_tempUnNE)))
                                                         +(1.0-phipar)*(piW(searcheffWW)*((1.0-piE(searcheffJ))*C_tempWWE+piE(searcheffJ)*((1.0-pLTpar)*C_tempWWE_JEi + pLTpar*C_tempWWE_JE))+
                                                                        (1.0-piW(searcheffWW))*(piE(searcheffJ)*((1.0-pLTpar)*C_tempUiE_JEi + pLTpar*C_tempUnE_JE) + (1.0-piE(searcheffJ))*((1.0-pLTpar)*C_tempUiE + pLTpar*C_tempUnE))));
#endif
#endif
        }
        if(typeU == 1) {    // not insured
            valfoc = utilc(cons,0) + disutilityW(searcheffWW) + disutilityE(searcheffJ) + betapar*(
                                                                                                   phipar*(piW(searcheffWW)*((1.0-piE(searcheffJ))*tempWWNE+piE(searcheffJ)*tempWWNE_JNE)+
                                                                                                           (1.0-piW(searcheffWW))*(piE(searcheffJ)*(tempUnNE_JNE) + (1.0-piE(searcheffJ))*(tempUnNE)))
                                                                                                   +(1.0-phipar)*(piW(searcheffWW)*((1.0-piE(searcheffJ))*tempWWE+piE(searcheffJ)*tempWWE_JE)+
                                                                                                                  (1.0-piW(searcheffWW))*(piE(searcheffJ)*(tempUnE_JE) + (1.0-piE(searcheffJ))*(tempUnE))));
#if OPT_compute_C_flow == 1
            focparamsU->C_val = utilc(cons,0) + betapar*(
                                                         phipar*(piW(searcheffWW)*((1.0-piE(searcheffJ))*C_tempWWNE+piE(searcheffJ)*C_tempWWNE_JNE)+
                                                                 (1.0-piW(searcheffWW))*(piE(searcheffJ)*(C_tempUnNE_JNE) + (1.0-piE(searcheffJ))*(C_tempUnNE)))
                                                         +(1.0-phipar)*(piW(searcheffWW)*((1.0-piE(searcheffJ))*C_tempWWE+piE(searcheffJ)*C_tempWWE_JE)+
                                                                        (1.0-piW(searcheffWW))*(piE(searcheffJ)*(C_tempUnE_JE) + (1.0-piE(searcheffJ))*(C_tempUnE))));
#endif
        }
    }
    
    valfoc=-valfoc;
    return valfoc;
    
}

